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A class of hybrid finite element methods 
for electromagnetics: a review 

John L. Volakis, A. Chatterjee and J. Gong 
Radiation Laboratory 

Department of Electrical Engineering and Computer Science 
The University of Michigan 
Ann Arbor MI 48109-2122 


Abstract 

Integral equation methods have generally been the workhorse for antenna and scat- 
tering computations. In the case of antennas, they continue to be the prominent 
computational approach, but for scattering applications the requirement for large- 
scale computations has turned researchers’ attention to near neighbor methods such 
as the finite element method, which has low O(N) storage requirements and is read- 
ily adaptable in modeling complex geometrical features and material inhomogeneities. 
In this paper, we review three hybrid finite element methods for simulating compos- 
ite scatterers, conformal microstrip antennas and finite periodic arrays. Specifically, 
we discuss the finite element method and its application to electromagnetic problems 
when combined with the boundary integral, absorbing boundary conditions and arti- 
ficial absorbers for terminating the mesh. Particular attention is given to large-scale 
simulations, methods and solvers for achieving low memory requirements and code 
performance on parallel computing architectures. 
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Introduction 


The 198(Ts brought to focus a variety of computational methods for electromagnetics applic- 
ations. Journal publications in this area grew dramatically and a variety of numerical codes 
were developed in industry and academia for analysis and design purposes. The methodo- 
logies used for developing these codes vary widely, but in general, they can be categorized 
in (a) moment method and (b) partial differential equation (PDE) formulations. Both were 
actually introduced in the sixties and seventies but the dramatic increase in computing re- 
sources over the past decade led to the development of accurate, useful and in many cases 
practical numerical codes for electromagnetic computations. Useful computer codes for wire 
(antennas and scatterers) and two-dimensional structures have been available since the sev- 
enties [1, 2, 3]. However, the more arduous task of developing useful codes applicable to 
complex three-dimensional structures took place, primarily, in the later part of the eighties 
[4]-[18] and continues at a rather fast pace. Driven by the capabilities of computing re- 
sources, a variety of basis functions, discretization schemes and solvers have been employed 
in connection with 3D moment method codes, all aimed at reducing the storage requirements 
and speeding-up the solution algorithms on either serial or parallel computers. 

To a great degree, though, leaps in computer speeds turned reseachers’ attention to 
method which have the lowest memory requirements, even at the expense of speed. PDE 
methods have inherently low memory requirements whether executed in time or frequency 
domain and have, thus, been courted by many research and engineering groups in industry, 
government and academia. Numerous successful applications of the Finite Difference-Time 
Domain (FD-TD) [1 1]— [13], [16], [19]— [21] and Finite Element Methods (FEM) [15, 17], 
[22]— [27] have been reported in the literature and by all accounts, the proliferation of the 
methods will continue to grow. Such methods have been routinely used for electrostatic and 
magnetostatic analyses, but during the past decade they were also developed for large scale 
computations involving hundreds of thousands of unknowns. 

Over the past few years, the group represented by these authors has been developing 
frequency domain PDE methods for scattering and antenna characterizations [17], [27]— [38] , 
In this paper, we review our progress with emphasis on three-dimensional applications. As 
is the case with any PDE method, the mesh truncation scheme plays an important role 
in assessing the accuracy and efficiency of the implementation. Various mesh termination 
schemes are mentioned here and their attributes are discussed. Generally, the mesh termin- 
ation scheme is either exact or approximate and each truncation scheme leads to a different 
hybrid FEM. Exact mesh truncation schemes have the drawback of destroying the sparsity 
of the finite element matrix but in most cases they also reduce the computational domain. 
Approximate truncation schemes rely on the use of absorbing boundary conditions or ma- 
terial absorbers whose purpose is to suppress wave reflections back into the computational 
domain. They retain the sparsity of the finite element matrix but one is required to increase 
the computational domain to ensure acceptable levels of accuracy. For antenna parameter 
prediction, there is a need for accurate near field calculations and because the absorbing 
boundary conditions (ABCs) lead to greater inaccuracies in the near zone, it is best that 
such predictions are carried out in connection with exact truncation schemes. 

Below, we first summarize the finite element formulation without reference to the mesh 
termination scheme. This formulation accounts for the presence of materials, conductors, 
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antenna feeds, resistive and impedance cards. The discretization of the weak form of the 
wave equation is then presented, followed by a discussion on the mesh truncation scheme and 
the properties of the resulting system. Both the boundary integral (BI) and ABC termina- 
tion schemes are presented with specific reference to the resulting system solvers. When a 
boundary integral is used to terminate the mesh, the iterative solver must be combined with 
the fast Fourier transform to maintain the storage requirements to 0(N). In the case of 
the finite element-ABC (FE-ABC) formulation, we discuss the code implementation and its 
execution and performance on parallel computing architectures. Finally, several examples 
are presented which demonstrate the utility and accuracy of the discussed hybrid FEM for- 
mulations for calculations pertaining to single element antennas, finite and infinite periodic 
arrays, and composite scatterers. Although this paper is certainly a review of hybrid finite 
element formulations, most of the included results have not appeared in the literature. 


Formulation 


The variational functional 


Consider the cavity configuration, shown in Figure 1, and the inhomogeneous scatterer shown 
in Figure 2. The cavity’s aperture shall be denoted by So whereas the same symbol will be 
used to denote the termination of the PDE’s computational domain in connection with the 
scatterer in Figure 2. As shown, the cavity is recessed in a ground plane and is intended to 
house the radiating elements and their feeding structure such as striplines, coaxial cables or 
microstrip lines. Any material filling can be handled by the presented formulation including 
resistive cards [34], lumped or distributed loads [33] often used for frequency and bandwidth 
control. Inhomogeneities are also assumed to be present in the dielectric scatterer of Figure 2. 

A goal with any PDE method, such as the finite element method, is the solution of the 
vector wave equation 


V x Q-V x Ej - ifcfcE = -jkoZoT - V x Q-M’j (1) 

in which E represents the total electric field, (e r ,^ r ) denote the relative permittivity and 
permeability of the computational domain, k 0 is the free space wavenumber and Zq is the 
free space intrinsic impedance. In addition, J* and M l denote the impressed electric and 
magnetic sources, respectively. 

Instead of solving (1) directly, for numerical accuracy and stability it is best to consider 
the weak form of (1). This is obtained by first multiplying (1) by a suitable weighting 
function and subsequently integrating over the computational domain [17]. An alternative 
approach is to consider the minimization of the functional [34] 


F(E) - 5 lllv f^ V * E ) ' ( V X E ) - *0CrE • E 


dv 


+ 


///„ 


E 


jk 0 Z 0 3 i + Vx — M' 
V/*r > 


dv 


3 


( 2 ) 


+ l -jk 0 Z 0 jj s ^(hxE)-(hxE)ds 
+ jk 0 Z 0 JJ E-(Hxn)ds 

in which R denotes the resistivity or impedance of the surface S r within So, and H is the 
corresponding magnetic field. We note that S r must be closed when it satisfies an opaque 
impedance boundary condition but can be open (i.e. a finite plate) if it satisfies the penetrable 
resistive sheet condition [39] 

hxnx E = -Rh x (H + - H) 

where H* denote the fields above and below the surface S r - As seen from (2), no special 
care is required for the evaluation of the integral over S r , in spite of the discontinuity in the 
magnetic field. The explicit knowledge of H is, however, required over the outer surface So 
for the construction of a discrete system of equations for the solution of E. Before proceeding 
with the latter, we note that for large computational domains it may be appropriate to re- 
express the functional (2) in terms of the scattered field. This reduces mesh propagation 
errors [40] and facilitates the implementation of the ABCs on So- Setting 


E = E' + E 1 


scat 


H = IT + H 


scat 


where E' denotes the excitation field impinging from the exterior of So and E scat is the 
scattered field, the variational functional F takes the form 


F(E ) 


- MIL 
III 


— (V x E scat ) • (V x E scat ) - klt r E 

fi r 


scat 


E 1 


scat 


dv 
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cat • \jkoZoT + V x j dv 

l -jk 0 Z 0 JJ s ±(n x E scat ) • (n x E scat ) ds 

jk 0 Z 0 jj s -^(n x E scat ) • (n x E ')ds 

JJJ v |-|-(V x E scat ) • (V x E*) - ^e r E scat • E' 

jk 0 Zo ff — E scat -(nxH *')ds 
JJs d fir 


dv 


+ \ ik ° Z °IJ s 

+ /( E') 


E 


scat 


So 


(h x H ) ds 


( 3 ) 


In this, /( E 1 ) is a function of the incident field whose specific form is of no consequence in the 
subsequent analysis. Most importantly, we observe in (3) the presence of the integrals over Sj 
and Vd , where Sd denotes the aggregate of all surfaces in V which coincide with a dielectric 
interface and Vd encompasses the volume of the enclosed dielectrics. The introduction of 
these integrals is required when working with the scattered field because at an interface 
E scat is associated with a jump in permeability. 
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Discrete system 

To construct a system of equations from (2), the volume V is first subdivided into a number 
of small elements such as those shown in Figure 3. Each of these small elements occupies a 
volume V e (e = 1, 2, 3, . . . , A/), where M denotes the total number of elements filling V. The 
field (total or scattered) in each element is then approximated by a linear or higher order 
expansion, viz. 

771 

E e = L = {W e } T {£ e } (4) 

J=1 

in which are the (vector) expansion basis functions [41]— [43] . With this representation, 
the unknown expansion coefficients E e 3 represent the field at the j th edge of the eth element 
having m sides/edges. The edge-based expansion (4) has zero divergence within each ele- 
ment as required by Maxwell’s equations and avoid an explicit specification of the field at 
metallic corners. Corresponding node-based field representations, which do not satisfy the 
zero divergence condition, may yield extraneous solutions unless augmented in some way 
by adding, for example, a penalty term to (2) or (3) [31], [44], or by instead working with 
potentials rather than the fields themselves [25]. Recently, however, node-based representa- 
tions have been proposed which remove this difficulty by enforcing the divergence condition 
in the process of constructing the node-based basis functions [45], 

The calculations to be presented here were carried out using a first order edge-based 
expansion in connection with tetrahedrals (also referred to as 1-form Witney elements [43]) 
[35, 37, 27], rectangular (right-angled) bricks [32, 33] or cylindrical-shell elements [46], Edge- 
based implementations using curvilinear bricks (also referred to as parabolic elements) have 
also been reported [47], [48]. It should be noted, nonetheless, that both edge-based and 
node-based expansion functions have their advantages and disadvantages [49]. Specifically, 
node-based functions are often found easier to implement and the first order element provides 
a linear representation of the field in all directions. However, apart from their difficulty with 
spurious solutions (a major restriction), they require double noding or other special care at 
dielectric interfaces and lead to larger matrix bandwidths. The edge-based elements do not 
require double-noding at dielectric interfaces, and usually lead to sparser matrices. However, 
the unknown count is generally doubled when using edge-based elements but this is balanced 
by the increased sparsity of the stiffness matrix [27], [49], Another drawback appears to be 
their inherent assumption that the field be constant along the edges forming the geometrical 
elements. 

The system of equations to be solved for E e - is obtained from (2) or (3) by a Rayleigh-Ritz 
procedure [22]. This amounts to differentiating the functional F with respect to each edge 
field and then setting it to zero. On substituting (4) into (2) or (3), taking the first variation 
of F with respect to E or E scat and assembling all elements, we obtain the system 

( BF" \ M M * N > m p 

^ = £[»']<«•} + £[£*]{£') + Els’) + Etc”} = o (s) 

\ UCj J e=l «=1 j=l p= 1 

The sum of the element submatrices [yl e ] results in a large sparse matrix and its elements 

are computed from the volume integral(s) over V in (2) or (3). The column { C p } results 
from the discretization of the volume integrals over Vj, the source volume V s and those 
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surface integrals over Sd and S r in (3) which involve the incident field. The matrix [5 s ] 
results from the discretization of the surface integral over S T into M„ elements, i.e. that 
over the impedance or resistive sheets. Finally, the column {<7*} gives the contribution of 
the boundary integral over So when the last is subdivided into N a surface elements. The 
computation of the matrix elements g s { depends on the representation of H (or H scat ) in 
terms of E (or E scat ) on So, and this will be addressed later. Specific expressions of the 
matrix elements appearing in (5) in connection with the functional (2) are 

A]j = JJJ V ' i(VxW').(VxW')-4W'.W' dv 

Cf = cf + c f 

+ jk 0 Z 0 // W (ff°xn) ds. 

J J s o 

and 

5?. = jkoZ 0 j f si i(n x W s ) • (n x W s ) ds 

9l = jk 0 Zo [[ W s • (H scat x n) ds. (8) 

'US* 

In these, H 5 ° represents the magnetic field on So when the antenna cavity is absent. Were 

we to consider the discretization of the functional in (3), the second integral in the definition 

of Cf in (7) must then be replaced by 

cf = jk 0 Z 0 // — W? • (n x H') ds + [f -t(n x W?) • (n x E*) ds 

+ /// p f— ( v x w ^) ■ ( v X E') - • E* dv. (9) 

JJJV* [/X r 

in which Vf is the volume of the pth element inside the dielectric. For antenna parameter 
computations, it is of course understood that (E*,H') are zero whereas for scattering com- 
putations (J',M') are zero. 

If the volume V contains lumped loads at the ith edge of the eth volume element, we 
must then add to the A\- matrix element the term 

Af - JJ E' W' • Wj dv (10) 

where l denotes the length of the edge load, Zi is the value of the load in Ohms, and s 
represents its cross-sectional area. Furthermore, for antenna parameter computations in 
connection with a coaxial cable feed, we must add to the left side of (5) the term [37] 

[D P ]{E P ] + {P p } 


( 6 ) 


( 7 ) 
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where 


(ii) 


and 



( 12 ) 


In these expressions, 5£ ab denotes the pth surface element which belong to S c& b as illustrated 
in Figure 1. The strength of the current in the center conductor of the cable is denoted by 7 0 
and p denotes the radial distance measured from this conductor. Further, p is the associated 
unit vector along the radial direction away from the center conductor. 

Below we consider the specifics associated with the computation of the column {# 5 }. It 
will be seen that different hybrid finite element methods result depending on the approach 
used for the computation of {<7 3 }. 


Finite element-boundary integral method 


The hybrid FE-BI method was introduced in the early seventies [50, 51] but its application to 
electromagnetic scattering did not occur until the eighties [52]— [55]. This method combines 
the adaptability of the FEM with the rigor of the BI or moment method. Because of 
the sparsity of the FEM matrix, this method is preferable to the volume moment method 
formulation [6], [7], which leads to completely full systems and is thus not attractive for 
iterative solution. The overall FE-BI system is partly sparse and partly full, but as will be 
seen later, the system can be efficiently solved with only O(N) memory demand on using 
uniform zoning on the mesh termination boundary which can be actually chosen to coincide 
with the target’s surface/aperture [17], [28]— [30] . 

In accordance with the FE-BI method, the magnetic field in the integral over So is 
represented by an integral equation which may take the form [17], [31], [32] 

H - H 9 ° + JJ s G(r, r') • [E(r') x z] ds' (13) 


where G is the dyadic Green’s function of the second kind such that n x (V x G) = 0 on Sq. 
For the antenna problem in Figure 1 , G becomes the free space dyadic Green’s function 


G = —j2k 0 Y 0 



e -J*o|r-r'| 

47r|r — r'| 


(14) 


with r and r' representing the observation and integration points, respectively, and I = 
xx + yy + zz is the unit dyad. The factor of 2 in (14) is due to image theory. In connection 
with this problem, i.e. that of a cavity recessed in a ground plane, H 5 ° is equal to the sum of 
the incident and reflected fields for scattering computations, or zero for antenna parameter 
computations. If So is not planar, then as noted earlier, for scattering computation H s ° is 
equal to the field scattered by an otherwise metallic surface So. That is, if So is cylindrical 
then H 9 ° is the total field due to a plane wave incident on a smooth metallic cylinder. 
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The integral in (13) can be discretized by introducing the expansion (for the s'th surface 
element) 

n 9 

(15) 


E(r') = £ Ef Wf(r') 


1=1 


where W* (r 7 ) are the basis for the volume element which borders the s'th surface element 
belonging to the surface So- As usual, n s denotes the number of edges forming the surface 
elements comprising Sq. When (15) is introduced in (13) and then into (8), we obtain 


N$ e n 3 

g‘ = EE EfGf/ 

S f = 1 1 = 1 


(16) 


where the coefficients Gfj S are, of course, dependent on the specific form of the Green’s 
function G(r, r 7 ). For an aperture recessed in a ground plane, G is given by (14), and on 
employing the divergence theorem twice, in this case we may explicitly write [31], [32] 


G''* = 2klJJ s \zxW‘ i (r)}yj S ',[wf( r ')xz}G 0 (Ty)dAds 


( 17 ) 


+ 


2 / Is- Vs ' [ W ^ (r) X *] ’ {/ L v# s ' [ W ’V) x *] G o(r, r') da ' j da 


when the normal to So is n = z. 

From (16), it is seen that the computation of the column {#} can be expressed in matrix 
form as 

{«'} = [ G]{E-'} (18) 

where the elements of the [ Q ] matrix are given by once the local indices of the last are 
transformed to global indices. Clearly, from (17), all elements of the [ Q ] matrix are non-zero, 
and in spite of the FE-BI method’s rigor, this destroys the sparsity of the overall system 
resulting from (5). Moreover, the memory requirement for storing [ Q] is O(iV^), where N se 
denotes the total number of edges used in the discretization of So. Unless special care is 
taken, a direct solution of (5) will be rather computationally intensive making the FE-BI 
method unsuitable for structures having electrically large So . As will be discussed later, if 
[Q] is made block Toeplitz, its memory demand is reduced to 0(N se ). Further, by using the 
fast Fourier transform (FFT) for computing {5 s }, the number of operations in carrying out 
the matrix product in (18) is also reduced from 0(N£ e ) to 0(N ae \og 2 N ae ) [57]. 


Finite Element-ABC Method 

In accordance with this method, some local boundary condition is enforced on the artificial 
surface So which encloses V . The goal with any ABC is to eliminate backward reflections 
from 5o, and a variety of these boundary conditions have been proposed over the years begin- 
ning with those of Bayliss and Turkel [58] and Enquist and Majda [59]. More recently, other 
ABCs such as those by Webb and Kanellopoulos [60] and Sesques [61] (see also Givoli [62]) 
have been proposed. All ABCs provide an approximate relation between the E and H fields 
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on the surface So which in most cases is derived by assuming a field expansion in inverse 
powers of r, where r is the radial distance from the center of Sq. If the ABC annihilates the 
first (2m + 1) inverse powers of r, it is then referred to as an mth order ABC. The zeroth 
order ABC is the Sommerfeld radiation condition 


- jk 0 Z 0 n x H scal = -jkoh xhx E acal (19) 

and a second order vector ABC proposed in [60] is given by 

-jk 0 Z 0 fi x H 8cal = oE 8Cat + /?V x [n (V x E 8cat ) J 

+ pv t (V • E 8cat ) (20) 

In this a = jk 0 , /? = l/(2jk 0 + 2/r), and the subscripts t and n denote the transverse and 
normal components to 5o, respectively. These ABCs were derived for a spherical surface, but 
have been found to work well when 5o is made piecewise planar so that it better conforms 
to the body’s surface [27]. In that case, the coefficient P reduces to /? = \/2jka. A second 
order ABC which allows for arbitrary surface curvatures was recently proposed [61], and it 
is given by 


- jk 0 Z 0 h x H = 


2kl E, - C • E 8 


iscat 


1 + JCm/ko 


2 scat 


2 (j k 0 - 2cm ) 2(jk 0 - 2 c m ) 

where c m denotes the mean curvature of the mesh truncation surface, 

gscat = v x j^ v x E scat) n J + y^y . E *cat) 


( 21 ) 


and 


l-UE?* = (3c2 +ci) 
2 1 2 


C2 + C 

-jko + - 1 - 1 


l ) 


U 


in which <^2 are the principal surface directions and are the corresponding surface 
curvatures. 

Regardless of which ABC is chosen, we shall denote the right hand sides of (19)— (21 ) by 
the function / 5 (E 8cat ), allowing us to rewrite (8) as 


9i — f f W- • P(E scat ) 

J J s* 


( 22 ) 


Then, upon introducing the expansion 15), we may express the column in matrix form 
as in (18). However, in this case all elements of [Q] vanish when s ^ s', and for s = s' we 
have that 

G” = // [aW^.W^^Vx w;) n .(VxW])„ 

- /3(V ■ W-,)(V • W*,)] is 

in connection with the ABC (20). As a result the [ Q ] matrix is very sparse, and this is the 
main advantage of the FE-ABC method. Obviously, the surface So must be placed at some 
distance away from the body, whereas in the case of the FE-BI method, So can be coincident 
with the outer surface of the structure. However, as will be demonstrated later, we have 
found that for scattering computations So can be placed only 0.3 wavelengths away from 
the target when the FE-ABC method is implemented in connection with the ABC given by 
( 20 ) [ 27 ], 
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Finite Element-Artificial Termination Structure 


Instead of enforcing an ABC on So, to eliminate backward reflections and ensure the outgoing 
nature of the waves, we may alternatively place in front of So some material absorber which 
serves the same purpose. This absorber may be of the same kind as that used in anechoic 
chambers, but to reduce its thickness to a minimum, it is instructive to design a thin multi- 
layer absorber having, possibly, artificial constitutive parameters. One such metal-backed 
mathematical absorber is illustrated in Figure 4 [36], 

The artificial termination structure (ATS) has the same advantages as the ABC in that 
it retains the sparsity of the overall system. Specifically, for metal-backed absorbers, the 
boundary integral over S 0 is zero, and thus the column {g 3 } in (5) is eliminated. So far, our 
implementations have showed that the ATS will not necessarily reduce the volume between 
So and the target. Moreover, the discretization of the artificial dielectric layers may cause 
some additional inconvenience. 


Mesh Generation 

A first step in any numerical analysis is the generation of an appropriate discrete model of the 
structure. For our case, this amounts to subdividing the volume enclosed by the surface So 
or the cavity into bricks or tetrahedrals, a procedure referred to as mesh generation. For 
rectangular cavities and structures, the mesh generation using rectangular bricks is fairly 
simple. An approach which was employed is to subdivide the volume V into layers of 
thickness Az. Each layer is then readily subdivided into N s x N s rectangular bricks. If the 
cavity housing the antenna structure is nonrectangular, rectangular bricks can still be used, 
but in that case staircasing will be unavoidable. However, for structures of arbitrary shape, 
it is necessary to resort to other elements such as tetrahedrals or curvilinear bricks. The 
generation of meshes using these elements requires sophisticated pre-processing packages 
and commercial software is available for this purpose. We have used a package marketed by 
SDRC I-DEAS [63] which employs a block meshing scheme and can be run on many different 
workstation platforms. A variety of other CAD packages are commercially available [64]. 

Regardless of which mesh generation scheme is employed, the final result is a file often 
called the Universal file. As a minimum, this file contains the following Tables: 

1. A list of the global volume element Nos. and their material code. 

2. A list of the nodes forming each volume element. 

3. A list of the (x,j/,z) coordinates for each node. 

4. A list of global edge Nos. and a list of the pairs of global node Nos. forming an edge. 

The last table is not usually contained in the Universal file but can be readily constructed 
from the first three tables. To the above lists we must also add information pertaining to 
the antenna geometry, feeding lines, resistive cards, lumped loads and about the presence 
of a coax cable or a probe feed. This information is generated by a preprocessor based on 
information provided by the user. For example, in specifying the location of a rectangular 
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patch or stripline, the user may enter the z location of the patch along with the ( x,y ) 
coordinates for two opposite corners. For non-planar metallic surfaces, it is necessary to flag 
or group the appropriate surface elements bordering the metallic surface. This may be done 
when generating the mesh. A similar procedure is required for identifying the presence of 
resistive cards. The resistivity must be specified for each resistive surface element, and this 
is conveniently done by referring to a materials identification code. 

The specification of lumped loads or conducting posts is relatively simple. These are 
assumed to be placed at the edges of the volume elements and consequently the mesh must 
be generated with this in mind. Once the mesh is completed, the user specifies the location 
or edge number and the impedance value of each load or post. If the post is perfectly 
conducting, i.e. a short, it is flagged as such and during matrix assembly, the field coefficient 
associated with the edge occupied by the short is set to zero. The center conductor of the 
coaxial cable is defined in this manner. 

Based on the structural and loading specifications, the following additional Tables may 
also be generated: 

• List of edge Nos. on a perfectly conducting surface (the fields of these edges are set to 
zero a priori when employing the total field formulation). 

• List of volume elements along with their global edge Nos. which border the outer 
surface. 

• List of volume elements and global edge Nos. bordering a resistive surface along with 
the associated resistivity. 

• List of global edge Nos. along with subglobal edge Nos. for combinations. 

The preprocessor can be run interactively, with or without a graphical user interface, de- 
pending on which mesh generation package is employed. The matrix elements are computed 
directly from the data provided by these tables using the algorithms outlined in the previous 
section. 


Matrix Assembly and System Solution 


Once all input tables are read, the goal is to assemble the matrix system 



[ 0 ] [ 0 ] 1 {#} _ {£} 
[0] [G]J {£(?} “ K) 


(23) 


where [.4] is generally a sparse square symmetric matrix, [Q\ was discussed earlier, the column 
{£T} denotes the field values at the edges not bordering So, and { E are the corresponding 
coefficients for the edges on So- The columns {6^} and {f>®} are computed from { Cf } 
as defined in (7) and (11). All matrix coefficients in the above system are numbered in 
accordance with the unique global or their own subglobal/local numbering schemes. If the 
pth global edge lies on a metallic surface or coincides with a metallic wire or post, E ^ or 
Ef? are set to zero or equal to the negative of the tangential incident field, depending on 
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which formulation is employed. Consequently the computation of A pn is skipped. Note 
that with this notation, m and n refer to global edge numbers, whereas the subscripts i 
and j used earlier are local indices for the eth element. Obviously, an edge may belong to 
more than one volume element and in that case A is constructed by summing the element 
equations for that edge. It is important to note that only the non-zero elements of the [A] 
matrix are stored, thus, ensuring the 0(N ) memory demand of the system. On the average, 
for tetrahedrals the minimum number of non-zero elements per row is 9 and the maximum 
number of nonzeros per row is 30. However, these numbers vary among geometries. In our 
implementations, the matrix was generated using a modified ITPACK storage scheme [65] 
and then stored as a long vector. This resulted in a storage requirement of approximately 
15AV to l6Nv, where Ny is the number of edges in V or unknowns. The ITPACK storage 
scheme is attractive for generating finite element matrices since the number of comparisons 
required while augmenting the matrix depends only on the locality of the corresponding edge 
and not on the number of unknowns. However, on using the ITPACK storage scheme for 
our application, almost half the space is lost in storing zeros. The modified ITPACK scheme 

[66] does alleviate this problem to a certain degree; however, 30% of the space is still lost in 
zero padding. The best trade-off between storage and speed for our application is obtained 
by storing the non-zero matrix elements in a long complex vector, the column indices in a 
long integer vector and the number of non-zeros per row in another integer vector. This data 
structure is referred to as the Compressed Sparse Row (CSR) format. In our implementation, 
a map of the number of non-zeros for each row is obtained through the modified ITPACK 
algorithm. The main program stores the matrix in CSR format, thus minimizing storage 
and sacrificing a bit of speed. The required storage is \§N complex words plus two integer 
vectors of length \hN and JV, respectively, which serve a s pointers. 

For the FE-ABC method, the [Q] matrix is also sparse and is generally combined with the 
[•4] matrix in solving (23). However, as noted earlier, [Q] is a full complex matrix when the 
FE-BI method is employed. In this case, if we were to solve (23) using a standard algorithm, 
it is appropriate to first rewrite it in the form 

{A v ){E:} + [Af}{E°} = {bl} (24) 

KHK} + [$}{£*} = {b»} (25) 

where [Af n ] and [A v ] are sparse matrices since they were decomposed from [.4], and gen- 
erally [Af] — [A®] 7 \ Also, the superscripts B and V imply fields/matrices associated with 

the boundary surface So or the volume V. The first of these can now be substituted into the 

second to obtain the system [52] 

{Mui -isiMfrv'']} {£„*■} = {* - [eiMfrMO ( 26 ) 

It is important to note that [Af ] is usually a symmetric, sparse matrix which is an important 
advantage in solving (26). We note that a substitution of (24) into (25) would require the 
inversion of a full complex matrix. Any of the many existing direct and iterative solvers 

[67] can be used for solving (26) including LU decomposition, banded matrix algorithm, 
and Lanczos algorithm. The last is better suited for large systems. Nevertheless, for a 
general full matrix, whichever solver is used will demand excessive computational intensity. 
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To alleviate this, it is important to make [Q] block Toeplitz by employing uniform gridding 
on the surface So- This is quite practical for flat surfaces (see Figure 1), and although it 
may cause some extra effort in creating the mesh and possibly staircasing, it is well worth 
it. Thus, if the volume element is a rectangular brick, then all elements on So should be 
chosen to be rectangular and equal. In the case of tetrahedrals, we must choose all surface 
elements to be right triangles of equal size. 

If [Q] is block Toeplitz, it can then be recast in a block circulant form [57], [68], [69], and 
in that case the FFT can be used to compute the product [G]{E‘ r f} in conjunction with an 
iterative solver. Specifically, if we denote a N tp x M tp circulant block of [C/] by [£? <p ], we may 
then write (m and n denote the x and y indices defining each element/edge) [57], [70]— [72] 

Mip-l N tp - 1 

E E Gr m ,, n _ n .<„ = FF T 2 -i{ FFT (GX,J F F T (£®j} (27) 

m'=0 n' = 0 

This operation requires that only one row of the [Q tp ] matrix block be stored reducing the 
storage of the overall [Q] matrix to 0(A r se ). When the biconjugate gradient (BCG) algorithm 
[73] for symmetric systems is used for the solution of (23), the overall storage requirement 
for rectangular elements is about ANy + 8jV se , excluding the storage needed for [A] and [Q]. 
For triangular surface elements the BCG algorithm’s storage requirement is about 4.5 times 
larger. The elements of [A] are extremely simple to compute and can be actually generated 
on the fly from the element submatrices for storage economy at the expense of computational 
efficiency. We have observed rather impressive convergence rates using this algorithm (with 
diagonal preconditioning) for FE-BI systems having more than 100,000 unknowns. For 
example, for 180,000 unknowns, convergence is achieved in 57 iterations. Another system 
of 25,000 unknowns required 66 iterations per angle with an average CPU time of 2 sec per 
iteration. Generally, the sampling rate was 12 to 15 elements per linear wavelength. 

The BCG algorithm was also used for solving the sparse system (23) in connection with 
the FE-ABC formulation. This formulation was exclusively used for treating targets asso- 
ciated with non-planar termination surfaces # 0 . Rather large-scale FE-ABC systems (over 
250,000 unknowns) were solved via the symmetric BCG algorithm with rather impressive 
convergence rates. Specifically, convergence was generally achieved with Nf 100 to Af/50 it- 
erations. The conjugate gradient squared (CGS) algorithm is generally faster than BCG but 
is more unstable since the residual polynomials are merely the squared BCG polynomials. 
Also, the CGS algorithm fails to exploit the symmetry of the matrix system. 

The performance of the BCG algorithm in solving the FE-ABC sparse system (23) was 
examined on the vector Cray YMP. It was observed that the vector updates and dot products 
reached speeds of about, 190 MFLOPS. However, the matrix vector products, which involve 
indirect addressing, ran at about 45.5 MFLOPS on a single processor of the 8-processor 
Cray YMP. As a rule of thumb, the BCG algorithm on the Cray YMP consumes about 
4.06 microseconds/iteration/unknown. A simple diagonal preconditioning can improve the 
convergence rate by about 35%, but a traditional ILU preconditioner was found less attractive 
because its convergence improvement was in most cases cancelled by the additional time 
required for its implementation. The parallelization of the BCG solver was also carried out 
on a 28-processor KSR1 parallel machine [66]. Figure 5 shows the speed-up for three versions 
of the sparse solver with diagonal preconditioning as the number of processors is increased 


13 


from 1 to 28. The maximum speed-up is observed to be about 19 and thus the solver was 
about 3.25 times faster on the 28-processor KSR1 than on the single-processor Cray YMP. 


Results 

In this section we present radar cross section (RCS) computations of cavity-backed antenna 
configurations, arrays (in a ground plane) and free-standing inhomogeneous structures. Once 
the fields in the computational domain are obtained, only the tangential fields are needed 
over a closed surface or the aperture So for an evaluation of the RCS or radiated fields. For 
free-standing structures, the far zone scattered field was obtained by applying the Stratton- 
Chu integral equation over the surface bounding the target. 

Example 1: RCS of a single patch residing in a rectangular cavity (FE-BI method) 

The geometry of this configuration is shown in Figure 6 along with the computed backscatter 
< Tee RCS pattern. The patch is 1.448" x 1.038" and resides on a dielectric substrate having 
thickness t = 0.057" and relative permittivity e r « 4.0. The substrate is housed in a 
2.89" x 2.10" x 0.057" rectangular cavity recessed in a ground plane. The excitation was a 6 - 
polarized plane wave at 9.2 GHz in the xz plane, and the analysis was done using rectangular 
bricks. As seen the computed (Tee pattern is in good agreement with the measured data [33]. 
We note that for this example the cavity terminations play a more important role because 
the computation frequency is not near the resonant frequency of the patch. The measured 
data were collected by placing the cavity-backed antenna configuration on a 5 foot long low 
cross section body [74]. 

Example 2: RCS of a circular patch near resonance (FE-BI method) 

The geometry of this patch configuration is shown in Figure 7 along with the calculated cr$e 
backscatter RCS as a function of frequency. The diameter of the circular patch is 1.3 cm 
and resides on a substrate 0.0787 cm thick, having e r « 2.33, The excitation is a plane 
wave incoming at (0{ = 60°, <f>i = 180°), and the computations were done using tetrahedral 
elements with the substrate terminated at a diameter of 1.6148 cm. However, because the 
calculations are near the patch resonance (~ 8 GHz), the substrate terminations do not 
contribute significantly to the RCS of the antenna. As seen, our computations are in good 
agreement with the measured data [75] which were collected with the patch placed in a cavity 
having a different periphery. 

Example 3: RCS of a patch with distributed resistive loading (FE-BI method) 

With the motivation of reducing the patch RCS at frequencies outside its operational band 
without substantial compromise in gain, we consider the placement of a narrow resistive 
ribbon (R-skirt) around the periphery of the patch as illustrated in Figure 8. Because the 
ribbon’s electrical width increases with increasing frequency, it is expected to have a rather 
noticeable effect on the RCS of the patch at higher frequencies. Thus, greater RCS reduction 
can be attained at higher frequencies without compromising the antenna’s gain performance. 
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If a certain level of RCS reduction is desired at resonance, lumped loads can still be used in 
conjunction with the R-skirt loading. 

The backscatter RCS data given in Figure 8 with and without loading correspond to a 
plane wave excitation incident in the xz plane and at an angle of 0 = 70° with respect to the 
z axis. However, as shown in Figure 9, the RCS reduction is fairly uniform over all 8 angles. 
As seen, the resistive skirt loading reduces the RCS in and out of the band as opposed to 
the narrowband performance of the lumped loads. Not surprisingly, the placement of the 
R-skirt is seen to slightly change the resonant frequency of the patch from 1.9522 GHz down 
to 1.744 GHz. From the inset in Figure 8, the RCS reduction at resonance is about 9 dB for 
all angles of incidences and thus the corresponding gain reduction is approximately 4.5 dB. 

Example 4: Scattering by an array of 3 x 30 slots (FE-BI method) 

Because of the low O(N) memory demand of the FE-BI method when combined with the 
FFT, one can consider the case of large finite arrays. The plot in Figure 10 is the 
backscatter RCS pattern of a 3 x 30 array of slots in a ground plane. Each slot has a 
rectangular 0.094" x 0.492" aperture and is 0.2059" deep. More than 100,000 unknowns 
were used for the simulation of this array and as seen the computed RCS pattern is in good 
agreement with corresponding moment method data. A characteristic of this pattern is the 
larger lobe at grazing which is often undesirable. As shown, this is completely suppressed 
when the last three slots from each end of the array are covered with a resistive card. The 
values of the resistive cards placed over the last three slots were 37.7 fi, 188.5 f l and 94.25 fi. 

Example 5: RCS of a composite cube (FE-ABC method) 

The geometry to be considered is an inhomogeneous cube 0.5A per side. Half of the cube 
(0.5A x 0.5A x 0.25A) consists of a dielectric having e r = 2 — j2 with its surfaces bordered by 
resistive cards having a resistivity equal to R = Zq. The other half of the cube is metallic. 
In Figure 11, we compare the computed principal plane a ^ RCS pattern with corresponding 
moment method data [76]. For the FE-ABC solution the scatterer was enclosed within a 
cubical outer boundary placed only 0.3A away from the scatterer. This resulted in a 30,000 
unknown system which converged in about 400 iterations on the Cray YMP when using 
the Sommerfeld radiation condition and in 1600 iterations when the second order ABC was 
used. For this geometry, the second order ABC did not provide a significant improvement 
in accuracy over the first order condition. 

Example 5: RCS of a 1.5A x 1A x 1A cavity (FE-ABC method) 

The parallelized version of the FE-ABC code was used to compute the RCS of 1.5A x 1A x 1A 
perfectly conducting inlet. In generating the mesh, this geometry was enclosed in a sphere 
of radius 1.35A. The resulting system had 224,476 unknowns and converged in about 5,000 
iterations on the KSR1. The computed backscatter pattern for the vertical polarization is 
shown in Figure 12. As seen, the computations are in good agreement with the measured 
data [77]. 
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Concluding Remarks 

The 1980’s brought to focus a variety of computational methods for electromagnetic ap- 
plications, including scattering, antennas, frequency selective surfaces, conformal arrays, 
microwave circuits and VLSI packaging. Not a single method is likely to predominate in 
all of these areas but it is nonetheless certain that future developments in computational 
electromagnetics will take into account the new architectures of high performance comput- 
ing facilities. The established integral equation methods will likely continue to be used for 
small and medium scale computations and optimized on high performance computing plat- 
forms. However, growing requirements to develop codes for modeling complex and large 
non-metallic geometrical structures; the availability of massively parallel computing archi- 
tectures and gallops in computer speeds do not seem to favor integral equation methods. 
For frequency domain computations, the finite element method, because of its geometrical 
adaptability, low O(N) memory demand and suitability for parallelization is undoubtedly 
a more promising computational method. The presented examples clearly demonstrated 
its accuracy for large scale electromagnetic computations but, needless to mention, there 
is much to be accomplished before the technique can be efficiently employed in modeling 
structures involving several millions of grid points. For example, there are needs for efficient 
generation algorithms, improved convergence solvers on parallel architectures, robust mesh 
truncation schemes, techniques for handling complex geometrical details and thin films, new 
adaptable element shapes and possibly new finite element formulations whose objective is 
to converge on a certain performance parameter (such as the RCS) rather than the fields 
themselves. 
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Figure 12. cr^e backscatter RCS pattern of a perfectly conducting rectangular in- 
let (1A x 1A x 1.5 A) taken in the yz plane. The black dots indicate 
data obtained from the KSR1 for the FE-ABC code and the solid line 
represents measured data [77]. 
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(a) Plan view 



Figure 1. Geometry of a cavity-backed antenna/array recessed in a ground plane. 





Figure 2. Illustration of a scatterer enclosed by an artificial surface 5 0 - 
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Figure 3. Geometries of various finite element shapes. 
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Figure 4. Scatterer enclosed by an artificial absorber termination, 
(a) Illustration of the enclose, (b) A specific ATS. 
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Speedup 



Figure 5. Speedup curves for the linear equation solver on the KSR1. 
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Figure 7. Comparison of measured and calculated backscatter RCS data for the 
shown circular patch as a function of frequency (0,* = 60°, <f>i = 180°). 
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Figure 8. <Tg$ backscatter RCS of the shown patch as a function of frequency with 
and without resistive skirt loading; incident wave direction = 70°, 
4>i = 180°. The patch is on a substrate of thickness 0.17558 cm having 
e r = 2.17. 
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Figure 9. backscatter RCS of the patch in Figure 8 as a function of angle 
with and without resistive loading. 
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Metal Boundary 



Figure 11. RCS pattern in the xz plane for the shown geometry. The solid curve 
is the FE-ABC pattern and the black dots are MoM data [76] for the 
E\ nc = 0 polarization. The specific geometry is a cube (a = b = 0.5A). 
Its lower half section is metallic and the upper half section is composed 
of dielectric (e r = 2 - j2, /i r = 1) bordered by a resistive card having 
R — Z 0 . 
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Figure 12. a$g backscatter RCS pattern of a perfectly conducting rectangular in- 
let (1A x 1A x L5A) taken in the yz plane. The black dots indicate 
data obtained from the KSR1 for the FE-ABC code and the solid line 
represents measured data [77]. 
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